/*
 * EBFV1_saturationTimeAdvance.cpp
 *
 *  Created on: 14/09/2012
 *      Author: rogsoares
 */

#include "EBFV1_hyperbolic.h"

namespace PRS{
	double EBFV1_hyperbolic::saturationTimeAdvance(TMesh* pTMesh_data, double delta_T){
		// take only nodes without wells
		int nrows, ncols, i, j, k;
		pTMesh_data->vertices_WithOutWells.getSize(nrows,ncols);
		for (int j=0; j<nrows; j++){
			VertexData* vdata = pTMesh_data->vertices_WithOutWells.getValue(j,0);
			vdata->Sw = vdata->Sw - (delta_T/vdata->weighted_volume)*vdata->nonvisc;				// calculate new saturation value

			// for visualization purposes only
			for (i=0; i<2; i++){
				vdata->vel[i] /= vdata->countedges;
			}
			vdata->countedges = 0;
			if ( vdata->Sw < 1.0e-6 ){
				vdata->Sw = .0;
			}
			if (vdata->Sw > 1.01 || vdata->Sw < .0){
				char text[256];
				sprintf(text,"Water saturation is out-of-bound [0 1]. Sw = %.10E on node %d\n",vdata->Sw,vdata->ID);
				throw Exception(__LINE__,__FILE__,text);
			}
		}

		// Calculate saturation advance for nodes with production wells
		double fw, Sw_old, wv, Sw_aux, Qwi, Sw0, dt_well, Qoi;
		if ( pSimPar->rankHasProductionWell() ){							// only ranks with production wells perform calculations below.
			const int N = pSimPar->getWellTimeDiscretion();												// number of fractional time-steps
			int numPW = pTMesh_data->getNumProductionWells();				// number of production wells
			for (int i=0; i<numPW; i++){
				double oilFlow = .0;
				pTMesh_data->productionWells_list[i].getSize(nrows,ncols);	// number of nodes for each production well
				for(int j=0; j<nrows; j++){
					VertexWell* vdata = pTMesh_data->productionWells_list[i].getValue(j,0);
					dt_well = (double)delta_T/N;						// time step for well nodes saturation
					Sw_old = vdata->vertex->Sw;
					wv = vdata->vertex->weighted_volume;
					Sw_aux = (dt_well/wv)*(vdata->vertex->nonvisc);
					// for visualization purposes only
					for (k=0; k<2; k++){
						vdata->vertex->vel[k] /= vdata->vertex->countedges;
					}
					vdata->vertex->countedges = 0;
					for (k=0; k<N; k++){
						Sw0 = Sw_old - Sw_aux;
						fw = pPPData->getFractionalFlux(Sw0);
						Qwi = fabs(fw*vdata->Qi);
						vdata->vertex->Sw = Sw0 - dt_well*(Qwi/wv);
						if ( vdata->vertex->Sw < 1.0e-6 ){
							vdata->vertex->Sw = .0;
						}
						Sw_old = vdata->vertex->Sw;
					}
					// calculate oil production
					Qoi = fabs(vdata->Qi) - fabs(fw*vdata->Qi);		// Oil production
					oilFlow += fabs(Qoi);
				}
				pOPManager->setRecoveredOilValue(oilFlow);
			}
		}

//		// print to file
//		ofstream fid;
//		fid.open("PRS.txt");
//		pTMesh_data->vertices_list[0].getSize(nrows,ncols);
//		for (int j=0; j<nrows; j++){
//			VertexData* vdata = pTMesh_data->vertices_list[0].getValue(j,0);
//			fid <<vdata->ID << "\t" << vdata->Sw << endl;
//		}
//		fid.close();
//		STOP();
		return 0;
	}
}
